* Michael Biggs, "The Dutch Protocol for Juvenile Transsexuals: Origins and Evidence",
* Journal of Sex and Marital Therapy, 2022; DOI 10.1080/0092623X.2022.2121238

* This Stata do-file compares Dutch and English adolescents suffering from gender dysphoria
* who were treated with GnRHa: 
* (1) 54 Dutch subjects, input manually from de Vries et al. 2011, table 2; 
* (2) 41 English subjects, downloaded from Carmichael et al. 2021, online supplement.

cd "/Users/michael/Documents/Gender/Health/Puberty blockers/Dutch & English 2"

**************************
* ENGLISH DATA 
**************************

use https://reshare.ukdataservice.ac.uk/854413/4/854413_GD_PublicDataset_v1.0.dta, clear
* From Carmichael, Polly, Gary Butler, Una Masic, Tim J. Cole, Bianca L. De Stavola, Sarah Davidson, Elin M. Skageberg, Sophie Khadr, and Russell Viner. 2021. `Short-Term Outcomes of Pubertal Suppression in a Selected Cohort of 12 to 15 Year Old Young People with Persistent Gender Dysphoria in the UK'. PLoS ONE 16:e0243894. doi: 10.1371/journal.pone.0243894.

foreach x in YSR CBCL { 
	foreach y in External Internal Total {
		gen `x'_`y'_end = `x'_12m_`y'_Tscore
		replace `x'_`y'_end = `x'_`y'_Tscore_24m if `x'_`y'_Tscore_24m!=.
		gen byte `x'_`y'_duration = cond(`x'_`y'_Tscore_24m!=., 24, cond(`x'_12m_`y'_Tscore!=., 12, .))
		gen `x'_`y'_change = `x'_`y'_end - `x'_b_`y'_Tscore
		ttest `x'_`y'_change==0
		}
	}
egen duration = rowmean(*_duration) 
tab duration
summ duration

di r(mean) / 12
summ *_change

save scratch, replace
local first = 1
foreach x in YSR CBCL { 
	foreach y in External Internal Total {
		di _n "`x'_`y'"
		use scratch, clear
		local collapse "(mean) t0_mean = `x'_b_`y'_Tscore (sd) t0_sd=`x'_b_`y'_Tscore"
		local collapse "`collapse' (mean) change_mean = `x'_`y'_change (sd) change_sd=`x'_`y'_change"
		local collapse "`collapse' (count) n = `x'_`y'_change"
		collapse `collapse'
		gen measure = "`x'_`y'"
		if `first' {
			save summary, replace
			local first = 0
			}
		else { 
			append using summary
			save summary, replace	
			}
		}
	}
	
gen sample = "English"
save summary, replace

**************************
* DUTCH DATA 
**************************

import delimited using rawdata_Dutch.csv, clear varnames(1) case(preserve)
* Data input manually from de Vries, Annelou L. C., Thomas D. Steensma, Theo A. H. Doreleijers, and Peggy T. Cohen‐Kettenis. 2011. `Puberty Suppression in Adolescents with Gender Identity Disorder: A Prospective Follow‐up Study'. Journal of Sexual Medicine 8:2276–83. doi: 10.1111/j.1743-6109.2010.01943.x.
append using summary
save summary, replace

**************************
* COMPARISON 
**************************

gen order = .
replace order = 6 if measure=="CBCL_Total"
replace order = 5 if measure=="CBCL_Internal"
replace order = 4 if measure=="CBCL_External"
replace order = 3 if measure=="YSR_Total"
replace order = 2 if measure=="YSR_Internal"
replace order = 1 if measure=="YSR_External"

**************

gen t0_lb = .
gen t0_ub = .
levelsof measure, local(mlevels)
foreach m of local mlevels {
	foreach c in Dutch English {
		foreach x in n t0_mean t0_sd  {
*			di "`x' `c' `m'"
			quietly summ `x' if sample=="`c'" & measure=="`m'"
			local `x'_`c' = r(mean)
			}
		quietly cii means `n_`c'' `t0_mean_`c'' `t0_sd_`c''
		quietly replace t0_lb = r(lb) if sample=="`c'" & measure=="`m'"
		quietly replace t0_ub = r(ub) if sample=="`c'" & measure=="`m'"			
		}
	di _n "`m'--`s':" 
	quietly sdtesti `n_Dutch' `t0_mean_Dutch' `t0_sd_Dutch' `n_English' `t0_mean_English' `t0_sd_English' 
	if r(p) < .05 {
		di "  (variance differs: p = " %6.5f (r(p))
		}
	local difference = 0
	quietly ttesti `n_Dutch' `t0_mean_Dutch' `t0_sd_Dutch' `n_English' `t0_mean_English' `t0_sd_English' `option'
	if r(p) < .05 {
		di "     p (equal)   = " %6.5f (r(p))
		local difference = `t0_mean_Dutch' - `t0_mean_English'
		}
	quietly ttesti `n_Dutch' `t0_mean_Dutch' `t0_sd_Dutch' `n_English' `t0_mean_English' `t0_sd_English', unequal
	if r(p) < .05 {
		di "     p (unequal) = " %6.5f (r(p)) 
		local difference = `t0_mean_Dutch' - `t0_mean_English'
		}
	if `difference' != 0	di "     Difference: " %5.3f (`difference')
	}	


**************************
* FIGURE 1: BASELINE 
**************************

local msize		"medlarge"
local fontsize  "12pt"
local Englishgrey "gs9"

gen row = order  //+ 6 + ((sex=="F") *8) + ((sample=="Dutch") *.25)
gsort +sample - row 
list measure order sample row 

forvalues r = 1/6 { 
	quietly levelsof measure_long if row==`r', local(label)
	local PFylabels "`PFylabels' `r' `label' " 
	}
di `"`PFylabels'"'

replace row = row + .15 if sample=="Dutch"
replace row = row - .15 if sample=="English"

twoway ///
	(rcap t0_lb t0_ub row if sample=="Dutch", horizontal lcolor(black) lwidth(medthick) msize(tiny)) ///
	(rcap t0_lb t0_ub row if sample=="English", horizontal color(`Englishgrey') lwidth(medthick) msize(tiny)) ///
	(scatter row t0_mean if sample=="Dutch", mcolor(black) msize(`msize')) ///
	(scatter row t0_mean if sample=="English", mcolor(`Englishgrey') msize(`msize')) ///
	, ///
	ytitle("") yscale(noline range(0.5)) ///
	ylabel( `PFylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("")   ///
	xlabel(50(10)70, grid labsize(`fontsize')) ///
	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l-4 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(3 4) label(3 "Amsterdam") label(4 "London") size(`fontsize')) ///
	name(Fig1, replace)
graph export Fig1.jpg, as(jpg) replace

******************************************************************************
******************************************************************************

**************************
* CHANGE OVER TIME
**************************

replace change_mean = (t1_mean - t0_mean) if sample=="Dutch" 

replace p = ".0005" if p=="<0.001"
destring p, replace

replace change_sd = abs(change_mean) * sqrt(n) / sqrt(F) if sample=="Dutch"
gen F_new = invFtail(1,n-1,p) if F<0
list measure F F_new if F<0
replace change_sd = abs(change_mean) * sqrt(n) / sqrt(F_new) if F<0

gen new_p = .
format p new_p %6.4f
gen change_lb = .
gen change_ub = .
local max = _N
forvalues i = 1/`max' {
	if change_sd[`i'] !=. {
		local n = n[`i']
		local change = change_mean[`i']
		local change_sd = change_sd[`i']
		quietly ttesti `n' `change' `change_sd' 0   
		quietly replace new_p = r(p) if _n==`i'
		quietly cii means `n' `change' `change_sd'
		quietly replace change_lb = r(lb) if _n==`i'
		quietly replace change_ub = r(ub) if _n==`i'			
		}
	}
list measure sample change_mean p new_p, clean noobs
format new_p %10.8f
list measure sample change_mean p new_p if new_p<0.0001, clean noobs

list measure t0_sd change_mean change_sd if F<0

**************************
* FIGURE 2: IMPROVEMENT
**************************

twoway ///
	(bar change_mean row if sample=="Dutch", horizontal color(gs3) barwidth(.25) ) ///
	(bar change_mean row if sample=="English", horizontal color(`Englishgrey') barwidth(.25)   ) /// 
	(rcap change_lb change_ub row if sample=="Dutch", horizontal lcolor(black) msize(tiny) lwidth(medthick)) /// 
	(rcap change_lb change_ub row if sample=="English", horizontal lcolor(`Englishgrey') msize(tiny) lwidth(medthick)) /// 
	, ///
	ytitle("") yscale(noline range(0.5)) xline(0, lcolor(black)) /// range(10.5 26)
	ylabel(`PFylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("", size(`fontsize'))   ///
	xlabel(, labsize(`fontsize') grid) ///
 	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l+10 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(1 2) label(1 "Amsterdam") label(2 "London") size(`fontsize')) ///
	name(Fig2, replace)
graph export Fig2.jpg, as(jpg) replace

foreach m of local mlevels {
	foreach c in Dutch English {
		foreach x in n change_mean change_sd  {
			quietly summ `x' if sample=="`c'" & measure=="`m'"
			local `x'_`c' = r(mean)
			}
		}
	di _n "`m'--`s':" 
	quietly sdtesti `n_Dutch' `change_mean_Dutch' `change_sd_Dutch' `n_English' `change_mean_English' `change_sd_English' 
	if r(p) < .05 {
		di "  (variance differs: p = " %6.5f (r(p))
		}
	local difference = 0
	quietly ttesti `n_Dutch' `change_mean_Dutch' `change_sd_Dutch' `n_English' `change_mean_English' `change_sd_English' `option'
	if r(p) < .05 {
		di "     p (equal)   = " %6.5f (r(p))
		local difference = `change_mean_Dutch' - `change_mean_English'
		}
	quietly ttesti `n_Dutch' `change_mean_Dutch' `change_sd_Dutch' `n_English' `change_mean_English' `change_sd_English', unequal
	if r(p) < .05 {
		di "     p (unequal) = " %6.5f (r(p)) 
		local difference = `change_mean_Dutch' - `change_mean_English'
		}
	if `difference' != 0	di "     Difference: " %5.3f (`difference')
	}	

